\documentclass[a4paper]{article}
\usepackage{xeCJK}
\setCJKmainfont{WenQuanYi Micro Hei}
\usepackage[affil-it]{authblk}
\usepackage[backend=bibtex,style=numeric]{biblatex}
\usepackage{graphicx}
\usepackage{geometry}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\geometry{margin=1.5cm, vmargin={0pt,1cm}}
\setlength{\topmargin}{-1cm}
\setlength{\paperheight}{29.7cm}
\setlength{\textheight}{25.3cm}

\addbibresource{citation.bib}

\begin{document}
% =================================================
\title{Numerical Analysis homework \# 6}

\author{王劼 Wang Jie 3220100105
  \thanks{Electronic address: \texttt{2645443470@qq.com}}}
\affil{(math), Zhejiang University }


\date{Due time: \today}

\maketitle

\begin{abstract}
    theoretical homework 6  
\end{abstract}





% ============================================
\section*{theoretical homework}

Complete the theoretical homework for section 5.6, and the programming homework is optional. Not submitting will not affect your regular grade, but submitting can be considered to offset the scores of previous homework assignments. \cite{wangheyu2024}

\subsection*{Question 5.6.1 I}

\textbf{Answer:}\\

\textbf{Part (a): Show that Simpson’s rule on \([-1, 1]\) with \(\rho = 1\) can be obtained by the expression}

Simpson’s rule for integration approximates the integral of a function \( y(t) \) by using a quadratic approximation (since \( p_3(t) \) is cubic) through the interval \([-1, 1]\). 
\[
\int_{-1}^{1} y(t) \, dt  = \int_{-1}^{1} (p_3(y; -1, 0, 0, 1; t) + \frac{y^{(4)}(\xi)}{4!}(t+1)t^2(t-1)) = \int_{-1}^{1} p_3(y; -1, 0, 0, 1; t) \, dt + E_S(y), \text{ while} \xi \in (-1,1).
\]
The error term, \( E_S(y) \), is the difference between the exact integral and the approximation given by the cubic interpolation. This is the key to the next part of the problem.

\textbf{Part (b): Derive \( E_S(y) \), the error in the approximation.}

To derive the error \( E_S(y) \), we use the fact from part(a):
\[
E_S(y) = \int_{-1}^{1} \frac{y^{(4)}(\xi)}{4!}(t+1)t^2(t-1) dt.
\]
The standard result for the error in a cubic interpolation is:
\[
E_S(y) = -\frac{1}{90} y^{(4)}(\xi),
\]
where \( y^{(4)}(\xi) \) is the fourth derivative of \( y(t) \).

\textbf{Part (c): Derive the composite Simpson’s rule and prove the theorem on its error estimation.}

The composite Simpson’s rule is based on dividing the interval \([a, b]\) into \( n \) subintervals (where \( n = 2m\) is even), and applying Simpson’s rule on each subinterval. This means we approximate each subinterval using a quadratic polynomial and sum up the results. The composite Simpson’s rule for an interval \([a, b]\) is given by:
\[
\int_a^b f(x) \, dx - S_n = \sum_{k=1}^m \left( \int_{x_{2k-2}}^{x_{2k}} f(x) \, dx - \frac{H}{6} [f(x_{2k-2}) + 4f(x_{2k-1}) + f(x_{2k})] \right) = -\frac{H^5}{2880} \sum_{k=1}^m f^{(4)}(\eta_k) \quad \text{where} \quad x_{2k-2} \leq \eta_k \leq x_{2k}
\]
By the continuity of \(f^{(4)}(x)\), there exists a point \(\eta \in [a, b]\) such that:
\[
\frac{1}{m} \sum_{k=1}^m f^{(4)}(\eta_k) = f^{(4)}(\eta).
\]
Thus, we obtain the error estimate for the composite Simpson's rule:
\[
R(f, S_n) = -\frac{mH^5}{2880} f^{(4)}(\eta) = -\frac{(b-a)}{2880} H^4 f^{(4)}(\eta) = - \frac{(b - a) h^4}{180} f^{(4)}(\eta), \quad a \leq \eta \leq b
\]
where \( H = 2h \).

\subsection*{Question 5.6.1 II}

\textbf{Answer:}\\

\subsection*{(a) Using the Composite Trapezoidal Rule}

The error for the composite trapezoidal rule is given by the formula:
\[
E_T = - \frac{(b - a) h^2}{12} f''(\xi),
\]
where \( \xi \in [a, b] \) and \( h \) is the width of each subinterval. The number of subintervals \( n \) is related to \( h \) by \( h = \frac{b - a}{n} \). We want the error to satisfy:
\[
|E_T| < 0.5 \times 10^{-6}.
\]
For the trapezoidal rule, the error bound is:
\[
|E_T| = \frac{(b - a) h^2}{12} \max_{x \in [0, 1]} |f''(x)|.
\]

\textbf{Step 1: Compute \( f''(x) \)}

For \( f(x) = e^{-x^2} \), we compute the first and second derivatives:
\[
f'(x) = -2x e^{-x^2}, \quad f''(x) = (4x^2 - 2) e^{-x^2}.
\]
The second derivative is maximized at \( x = 0 \), where:
\[
f''(0) = -2.
\]
Thus, \( \max_{x \in [0, 1]} |f''(x)| = 2 \).

\textbf{Step 2: Apply the Error Bound}

Substituting the values into the error bound, we get:
\[
|E_T| = \frac{h^2}{6}.
\]
We want this to be less than \( 0.5 \times 10^{-6} \), so:
\[
\frac{h^2}{6} < 0.5 \times 10^{-6}.
\]
\[
n \geq 578.
\]
Thus, the number of subintervals required by the composite trapezoidal rule is approximately \( n = 578 \).

\subsection*{(b) Using the Composite Simpson’s Rule}

The error for the composite Simpson’s rule is given by the formula:
\[
E_S = - \frac{(b - a) h^4}{180} f^{(4)}(\xi),
\]
where \( \xi \in [a, b] \) and \( h = \frac{b - a}{n} \). We want the error to satisfy:
\[
|E_S| < 0.5 \times 10^{-6}.
\]

\textbf{Step 1: Compute \( f^{(4)}(x) \)}

For \( f(x) = e^{-x^2} \), we compute the first few derivatives:
\[
f'(x) = -2x e^{-x^2}, \quad f''(x) = (4x^2 - 2) e^{-x^2}, \quad f^{(3)}(x) = (12x - 8x^3) e^{-x^2},
\]
\[
f^{(4)}(x) = (16x^4 - 48x^2 + 12) e^{-x^2}.
\]
The fourth derivative is maximized at \( x = 0 \), where:
\[
f^{(4)}(0) = 12.
\]
Thus, \( \max_{x \in [0, 1]} |f^{(4)}(x)| = 12 \).

\textbf{Step 2: Apply the Error Bound}

Substituting the values into the error bound, we get:
\[
|E_S| = \frac{h^4}{15}.
\]
We want this to be less than \( 0.5 \times 10^{-6} \), so:
\[
\frac{h^4}{15} < 0.5 \times 10^{-6}.
\]
\[
n \geq 20.
\]
Thus, the number of subintervals required by the composite Simpson’s rule is approximately \( n = 20 \).

\section*{Summary}

\begin{itemize}
  \item (a) The number of subintervals required by the composite trapezoidal rule is approximately \( n = 578 \).
  \item (b) The number of subintervals required by the composite Simpson’s rule is approximately \( n = 20 \).
\end{itemize}

\subsection*{Question 5.6.1 III}

\textbf{Answer:}\\

\subsection*{(a) Construct a polynomial \(\pi_2(t) = t^2 + at + b\) that is orthogonal to \(P_1(t) = t\) with respect to the weight function \(\rho(t) = e^{-t}\)}

We are given that \(\pi_2(t)\) is a quadratic polynomial, so assume the form:

\[
\pi_2(t) = t^2 + at + b.
\]

We need to find the constants \(a\) and \(b\) such that:

\[
\int_0^\infty t \cdot (t^2 + at + b) e^{-t} \, dt = 0.
\]

\[
\int_0^\infty (t^2 + at + b) e^{-t} \, dt = 0.
\]

So we will solve by \( \int_0^\infty t^m \cdot e^{-t} \, dt = m!\):

\begin{align*}
  6 + 2a + b &= 0 \\
  2 + a + b &= 0
\end{align*}

We have:

\[
a=-4, b=2.
\]

Thus, the polynomial \(\pi_2(t)\) is:

\[
\pi_2(t) = t^2 - 4t + 2.
\]

\subsection*{(b) Derive the two-point Gauss-Laguerre quadrature formula}

\[
\int_0^\infty f(t) e^{-t} \, dt = \int_0^\infty H_3(t) e^{-t} \, dt + \int_0^\infty \frac{f^{(4)}(\tau)}{4!} w^2(t) e^{-t} \, dt = w_1 f(t_1) + w_2 f(t_2) + E_2(f),
\]

The error term \(E_2(f)\) for Gauss-Laguerre quadrature is expressed as:

\[
E_2(f) = \int_0^\infty \frac{f^{(4)}(\tau)}{4!} w^2(t) e^{-t} \, dt = \frac{f^{(4)}(\tau)}{6}.
\]
where \(f^{(4)}(\tau)\) is the fourth derivative of \(f(t)\).

\subsection*{(c) Apply the formula to approximate \(I = \int_0^\infty \frac{1}{1 + t} e^{-t} \, dt\)}

We compute:

\[
t_1 = 2 - \sqrt{2} \approx 0.5857864, t_2 = 2 + \sqrt{2} \approx 3.4142136.
\]

We konw \(w_1 \approx 0.8535534, w_2 \approx 0.1464466\), so we compute:

\[
w_1 f(t_1) + w_2 f(t_2) \approx 0.5714285877181063
\]

The error term is:

\[
E_2(f) = \frac{1}{(1+\tau)^5} \int_0^\infty (t^2-4t+2)^2`' e^{-t} \, dt = \frac{4}{(1+\tau)^5} \approx 0.024918773281893625
\]

Thus, the estimated error is about \(0.024918773281893625\), and the total approximation for the unknown quantity \(\tau\) is:

\[
\tau \approx 1.757807000787479.
\]

\subsection*{Question 5.6.1 IV}

\textbf{Answer:}\\

\section*{(a) Find \( h_m(t) \) and \( q_m(t) \) in the form}
We know:
\[
h_m(x_m)=1,
q_m(x_m)=0,
h_m'(x_m)=0,
q_m'(x_m)=1.
\]

which equals to:
\begin{align*}
a_m+b_mx_m &= 1 \\
c_m+d_mx_m &= 0 \\
b_m+2(a_m+b_mx_m)\sum_{i \neq m}^{n} \frac{1}{x_m-x_i} &= 0 \\
c_m+2(c_m+d_mx_m)\sum_{i \neq m}^{n} \frac{1}{x_m-x_i} &= 1 
\end{align*}

So we get:
\[
a_m=1-\sum_{i \neq m}^{n} \frac{2x_m}{x_i-x_m},
b_m=\sum_{i \neq m}^{n} \frac{2}{x_i-x_m},
c_m=-x_m,
d_m=1,
\]

\section*{(b) Deriving the Quadrature Rule}

Obviously \(p(t) = \sum_{m=1}^n \left[ h_m(t) f_m + q_m(t) f'_m \right] \in P_{2n-1}\) is Hermite interpolation polynomial on \(x_1, x_2, \dots, x_n\).

So we get:
\[
I_n(f) = \sum_{k=1}^n \left[ \int_a^b h_k(t) \rho(t) \, dt f(x_k) + \int_a^b q_k(t) \rho(t) \, dt f'(x_k) \right],
\]

which equals to\(w_k=\int_a^b \left[1-\sum_{i \neq k}^{n} \frac{2x_k}{x_i-x_k} + \sum_{i \neq k}^{n} \frac{2}{x_i-x_k} t \right]l_k^2(t) \rho(t) \, dt, \mu_k=\int_a^b (t-x_k)l_k^2(t) \rho(t) \, dt\)

It is exact for \(p \in P_{2n-1}\), so it satisfies \(E_n(p) = 0\).

\section*{(c) Conditions for \( \mu_k = 0 \)}

From (a) and (b), we know \(\mu_k=\int_a^b q_k(t) \rho(t) \, dt=\int_a^b (t-x_k)l_k^2(t) \rho(t) \, dt\)

When \(\mu_k = 0\) we get:
\[
t-x_k \text{ is orthogonal to } l_k^2(t),
\]
which \(k = 1, 2, \cdots, n\)


% ===============================================
\section*{ \center{\normalsize {Acknowledgement}} }
None.


\printbibliography

\end{document}